Questions:
Objectives:
dplyr and tidyr packages.Keypoints:
dplyr and tidyrBracket subsetting is handy, but it can be cumbersome and difficult to read, especially for complicated operations.
Some packages can greatly facilitate our task when we manipulate data.
Packages in R are basically sets of additional functions that let you
do more stuff. The functions we’ve been using so far, like str() or
data.frame(), come built into R; Loading packages can give you access to other
specific functions. Before you use a package for the first time you need to install
it on your machine, and then you should import it in every subsequent
R session when you need it.
The package dplyr provides powerful tools for data manipulation tasks.
It is built to work directly with data frames, with many manipulation tasks
optimised.
As we will see latter on, sometimes we want a data frame to be reshaped to be able
to do some specific analyses or for visualisation. The package tidyr addresses
this common problem of reshaping data and provides tools for manipulating
data in a tidy way.
To learn more about dplyr and tidyr after the workshop,
you may want to check out this handy data transformation with
dplyr
cheatsheet
and this one about
tidyr.
tidyverse package is an “umbrella-package” that installs
several useful packages for data analysis which work well together,
such as tidyr, dplyr, ggplot2, tibble, etc.
These packages help us to work and interact with the data.
They allow us to do many things with your data, such as subsetting, transforming,
visualising, etc.To install and load the tidyverse package type:
BiocManager::install("tidyverse")
# load the tidyverse packages, incl. dplyr
library("tidyverse")
rna <- read_csv("../data/rnaseq.csv")
# view the data
rna
Notice that the class of the data is now referred to as a “tibble”.
Tibbles tweak some of the behaviors of the data frame objects we introduced in the previously. The data structure is very similar to a data frame. For our purposes the only differences are that:
It displays the data type of each column under its name.
Note that <dbl> is a data type defined to hold numeric values with
decimal points.
It only prints the first few rows of data and only as many columns as fit on one screen.
We are now going to learn some of the most common dplyr functions:
select(): subset columnsfilter(): subset rows on conditionsmutate(): create new columns by using information from other columnsgroup_by() and summarise(): create summary statistics on grouped dataarrange(): sort resultscount(): count discrete valuesTo select columns of a data frame, use select(). The first argument
to this function is the data frame (rna), and the subsequent
arguments are the columns to keep.
select(rna, c(gene, sample, tissue, expression))
# select(rna, c("gene", "sample", "tissue", "expression"))
To select all columns except certain ones, put a “-” in front of the variable to exclude it.
select(rna, -tissue, -organism)
This will select all the variables in rna except tissue
and organism.
To choose rows based on a specific criteria, use filter():
filter(rna, sex == "Male")
filter(rna, sex == "Male" & infection == "NonInfected")
Now let’s imagine we are interested in the human homologs of the mouse
genes analysed in this dataset. This information can be found in the
last column of the rna tibble, named
hsapiens_homolog_associated_gene_name. To visualise it easily, we
will create a new table containing just the 2 columns gene and
hsapiens_homolog_associated_gene_name.
genes <- select(rna, gene, hsapiens_homolog_associated_gene_name)
genes
Some mouse genes have no human homologs. These can be retrieved using
filter() and the is.na() function, that determines whether
something is an NA.
filter(genes, is.na(hsapiens_homolog_associated_gene_name))
If we want to keep only mouse genes that have a human homolog, we can
insert a “!” symbol that negates the result, so we’re asking for
every row where hsapiens_homolog_associated_gene_name is not an
NA.
filter(genes, !is.na(hsapiens_homolog_associated_gene_name))
What if you want to select and filter at the same time? One way to accomplish this is using pipes, which are a recent addition to R. Pipes let you take the output of one function and send it directly to the next, which is useful when you need to do many things to the same dataset.
Pipes in R look like %>%.
In the above code, we use the pipe to send the rna dataset first
through filter() to keep rows where sex is Male, then through
select() to keep only the gene, sample, tissue, and
expressioncolumns.
The pipe %>% takes the object on its left and passes it directly as
the first argument to the function on its right, we don’t need to
explicitly include the data frame as an argument to the filter() and
select() functions any more.
rna %>%
filter(sex == "Male") %>%
select(gene, sample, tissue, expression)
Some may find it helpful to read the pipe like the word “then”. For instance,
in the above example, we took the data frame rna, then we filtered
for rows with sex == "Male", then we selected columns gene, sample,
tissue, and expression.
The dplyr functions by themselves are somewhat simple, but by
combining them into linear workflows with the pipe, we can accomplish
more complex manipulations of data frames.
If we want to create a new object with this smaller version of the data, we can assign it a new name:
rna3 <- rna %>%
filter(sex == "Male") %>%
select(gene, sample, tissue, expression)
rna3
Using pipes, subset the rna data to keep observations in female mice at time 0,
where the gene has an expression higher than 50000, and retain only the columns
gene, sample, time, expression and age.
Frequently you’ll want to create new columns based on the values of existing
columns, for example to do unit conversions, or to find the ratio of values in two
columns. For this we’ll use mutate().
To create a new column of time in hours:
rna %>%
mutate(time_hours = time * 24) %>%
select(time, time_hours)
You can also create a second new column based on the first new column within the same call of mutate():
rna %>%
mutate(time_hours = time * 24,
time_mn = time_hours * 60) %>%
select(time, time_hours, time_mn)
Many data analysis tasks can be approached using the
split-apply-combine paradigm: split the data into groups, apply some
analysis to each group, and then combine the results. dplyr
makes this very easy through the use of the group_by() function.
rna %>%
group_by(gene)
The group_by() function doesn’t perform any data processing, it
groups the data into subsets: in the example above, our initial
tibble of 32428 observations is split into
1474 groups based on the gene variable.
We could similarly decide to group the tibble by the samples:
rna %>%
group_by(sample)
Here our initial tibble of 32428 observations is split into
22 groups based on the sample variable.
Once the data has been grouped, subsequent operations will be applied on each group independently.
summarise() functiongroup_by() is often used together with summarise(), which
collapses each group into a single-row summary of that group.
group_by() takes as arguments the column names that contain the
categorical variables for which you want to calculate the summary
statistics. So to compute the mean expression by gene:
rna %>%
group_by(gene) %>%
summarise(mean_expression = mean(expression))
We could also want to calculate the mean expression levels of all genes in each sample:
rna %>%
group_by(sample) %>%
summarise(mean_expression = mean(expression))
But we can can also group by multiple columns:
rna %>%
group_by(gene, infection, time) %>%
summarise(mean_expression = mean(expression))
## `summarise()` has grouped output by 'gene', 'infection'. You can override using
## the `.groups` argument.
Once the data is grouped, you can also summarise multiple variables at the same
time (and not necessarily on the same variable). For instance, we could add a
column indicating the median expression by gene and by condition:
rna %>%
group_by(gene, infection, time) %>%
summarise(mean_expression = mean(expression),
median_expression = median(expression))
## `summarise()` has grouped output by 'gene', 'infection'. You can override using
## the `.groups` argument.
Calculate the mean expression level of gene “Dok3” by timepoints.
In the rna tibble, the rows contain expression values (the unit) that are
associated with a combination of 2 other variables: gene and sample.
All the other columns correspond to variables describing either the sample (organism, age, sex, …) or the gene (gene_biotype, ENTREZ_ID, product, …). The variables that don’t change with genes or with samples will have the same value in all the rows.
rna %>%
arrange(gene) #sorts the dataframe
This structure is called a long-format, as one column contains all the values,
and other column(s) list(s) the context of the value.
In certain cases, the long-format is not really “human-readable”, and another format,
a wide-format is preferred, as a more compact way of representing the data.
This is typically the case with gene expression values that scientists are used to
look as matrices, were rows represent genes and columns represent samples.
In this format, it would therefore become straightforward to explore the relationship between the gene expression levels within, and between, the samples.
To convert the gene expression values from rna into a wide-format,
we need to create a new table where the values of the sample column would
become the names of column variables.
The key point here is that we are still following a tidy data structure, but we have reshaped the data according to the observations of interest: expression levels per gene instead of recording them per gene and per sample.
The opposite transformation would be to transform column names into values of a new variable.
We can do both these of transformations with two tidyr functions,
pivot_longer() and pivot_wider() (see
here for
details).
Let’s select the first 3 columns of rna and use pivot_wider()
to transform the data into a wide-format.
rna_exp <- rna %>%
select(gene, sample, expression)
rna_exp
pivot_wider takes three main arguments:
names_from : the column whose values will become new column
names;values_from: the column whose values will fill the new
columns.
Figure 1: Wide pivot of the rna data
rna_wide <- rna_exp %>%
pivot_wider(names_from = sample,
values_from = expression)
rna_wide
Note that by default, the pivot_wider() function will add NA for missing values.
Let’s imagine that for some reason, we had some missing expression values for some genes in certain samples. In the following fictive example, the gene Cyp2d22 has only one expression value, in GSM2545338 sample.
rna_with_missing_values
By default, the pivot_wider() function will add NA for missing
values. This can be parameterised with the values_fill argument of
the pivot_wider() function.
rna_with_missing_values %>%
pivot_wider(names_from = sample,
values_from = expression)
rna_with_missing_values %>%
pivot_wider(names_from = sample,
values_from = expression,
values_fill = 0)
In the opposite situation we are using the column names and turning them into a pair of new variables. One variable represents the column names as values, and the other variable contains the values previously associated with the column names.
pivot_longer() takes four main arguments:
names_to: the new column name we wish to create and populate with the
current column names;values_to: the new column name we wish to create and populate with
current values;names_to and
values_to variables (or to drop).
Figure 2: Long pivot of the rna data
To recreate rna_long from rna_wide we would create a key
called sample and value called expression and use all columns
except gene for the key variable. Here we drop gene column
with a minus sign.
Notice how the new variable names are to be quoted here.
rna_long <- rna_wide %>%
pivot_longer(names_to = "sample",
values_to = "expression",
-gene)
rna_long
Note that if we had missing values in the wide-format, the NA would be
included in the new long format.
Remember our previous fictive tibble containing missing values:
rna_with_missing_values
wide_with_NA <- rna_with_missing_values %>%
pivot_wider(names_from = sample,
values_from = expression)
wide_with_NA
wide_with_NA %>%
pivot_longer(names_to = "sample",
values_to = "expression",
-gene)
Now that you have learned how to use dplyr to extract information from
or summarise your raw data, you may want to export these new data sets to share
them with your collaborators or for archival.
Similar to the read_csv() function used for reading CSV files into R, there is
a write_csv() function that generates CSV files from data frames.
Before using write_csv(), we are going to create a new folder, data_output,
in our working directory that will store this generated dataset. We don’t want
to write generated datasets in the same directory as our raw data.
It’s good practice to keep them separate. The data folder should only contain
the raw, unaltered data, and should be left alone to make sure we don’t delete
or modify it. In contrast, our script will generate the contents of the data_output
directory, so even if the files it contains are deleted, we can always
re-generate them.
Let’s use write_csv() to save the rna_wide table that we have created previously.
write_csv(rna_wide, file = "data_output/rna_wide.csv")